High frequency sound waves in vitreous silica 
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We report a molecular dynamics simulation study of the sound waves in vitreous silica in the meso- 
scopic exchanged momentum range. The calculated dynamical structure factors are in quantitative 
agreement with recent experimental inelastic neutron and x-ray scattering data. The analysis of the 
longitudinal and transverse current spectra allows to discriminate between opposite interpretations 
of the existing experimental data in favour of the propagating nature of the high frequency sound 
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PACS numbers : 61.43. Bn, 61.43.Fs, 63.50.+X 

In the last 20 years great effort has been devoted to 
understanding the microscopic dynamics in topologically 
disordered systems. The need originates from the expe- 
rimental observation that a large variety of glasses ex- 
hibits common dynamical and thermodynamical proper- 
ties, which are anomalous with respect to the correspon- 
ding crystals Q. In particular, around 5 K the ther- 
mal conductivity has a plateau and in the 10 4 30 K 
temperature range the reduced specific heat Cp(T)/T 3 
goes through a maximum, indicating the existence of ex- 
cess vibrations in the low-frequency range with respect 
to the Debye model. In a more extended range, below 
Rj 10 meV, inelastic incoherent neutron Q and Raman 
j| scattering reveal a strong broad band, usually called 
boson peak, which can again be related to an excess of the 
vibrational density of states. At the present, a unique mi- 
croscopic description of these excess modes does not ex- 
ist and their theoretical interpretations include localized 
vibrational states B, soft anharmonic vibrations re- 
laxational motions [Q] , intrinsically diffusive || and prop- 
agating modes. 

Among disordered systems, vitreous silica is a typical 
strong glass former and has widely been studied, both ex- 
perimentally |^,^,^;D and numerically ||. The temper- 
ature analysis of the inelastic neutron-scattering inten- 
sity H suggested the coexistence of sound waves, below 
« 4 meV, with a second class of local harmonic exci- 
tations, which become increasingly anharmonic towards 
the lowest frequencies. However, kinematic reasons did 
not (and do not yet) allow the use of the neutron scat- 
tering technique to directly study the acoustic excita- 
tions in the boson peak energy region. Only recently it 
became possible to utilize the inelastic x-ray scattering 
(IXS) with meV energy resolution [[ll] to measure the 
dynamic structure factor, S(Q,E), in the first pseudo- 
Brillouin zone. Using this technique, the high frequency 
dynamics of v-Si02 at T = 1050 K has been studied in 
the 14-6 nm -1 momentum transfer range. The measured 
S(Q, E) clearly indicate the existence of collective propa- 
gating excitations up to Q = 4 nm -1 , i.e. up to an energy 
of ps 13 meV. The energy of these excitations extends 



both below and above the boson peak region, located at 
£«4 meV [Q . Recently this result has been confirmed 
in the whole 1800 4 300 K temperatures range p2| . At 
variance with this interpretation, analogous room tem- 
perature IXS data combined with non-kinematic inelastic 
neutron (INS) spectra Jl3| were interpreted as an indi- 
cation of the existence of a cross-over from sound waves 
to localized acoustic modes, which should be responsible 
for the boson peak. 

In order to clarify the propagating nature of low fre- 
quency modes around the boson peak region, this letter 
presents a molecular dynamics (MD) and normal mode 
analysis (NMA) study of the high frequency dynamics in 
vitreous Si02 in the Q — 1.4 4 5.1 nm -1 range. The 
reliability of our investigation is validated by the good 
agreement of the calculated S(Q,E) with the available 
INS § and IXS § data. The detailed knowledge of the 
spectral shape allows us to discriminate between the two 
proposed interpretations (?],[lj| and to establish the prop- 
agating nature of both longitudinal and transverse modes 
giving rise to the boson peak. Moreover, the widths of 
the excitations calculated in the harmonic approxima- 
tion (i.e. via NMA) are in quantitative agreement with 
both MD and experimental data, indicating that in the 
examined Q range their origin is structural rather than 
dynamical. 

Standard microcanonical MD simulations with peri- 
odic boundary conditions were carried out for systems 
consisting of 648 and 5184 particles enclosed in cubic 
boxes of length L = 2.139 nm and L = 4.278 nm re- 
spectively (mass density p = 2.20 g/cm 3 ). We used the 
two- and three-body interaction potential proposed by 
Vashista et al. [|l4| which, compared with the simpler- 
to-deal two-body potentials, gives a better description of 
the structural properties of Si02. Details of this study 
will be discussed in a subsequent paper. The electro- 
static long range interactions were taken into account by 
the tapered reaction field method and the equations of 
motion were integrated with leap-frog algorithm using a 
time step of At — 0.5 fs. A /3-cristobalite crystal was 
melted at 15, 000 K and equilibrated for a long time 
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so that the initial state has no effect on the final con- 
figuration. The system was then cooled and thermal- 
ized through different temperatures ranging from 5000 
K down to 300 K. The temperature dependence of the 
total internal energy and of the self-diffusion constants 
Dsi, Do are studied to recognise the structural arrest, 
which is found at about 2700 K, a value definitely higher 
than the experimental glass-transition temperature |l5| ] 
(T g ps 1450 K). This can be an effect, at least partially 
of the quench rate [pi, which in our MD simulation is 
much higher (« lO^^^C/s) than the experimental value 
(w 10 2 K/s). Two different techniques were adopted to 
study the dynamic of the system in the glassy state: i) 
(MD) the finite temperature (T = 1500, 1000, 600, and 
300 K) atomic dynamics has been followed for 60 ps by 
recording the instantaneous configurations every 0.01 ps; 
ii) (NMA) the dynamical matrix has been calculated in 
a " glassy" (metastable) energy minimum, reached by the 
steepest-descent method [[l7|, and then either diagonal- 
ized to obtain a complete set of eigenvalues and eigenvec- 
tors or used in the application of the spectral moments 
method ||,||. 

We investigated the dynamical correlation in v-SiC>2 
in terms of neutron- (Sn(Q, E)) and x-ray- (Sx{Q, E)) 
weighted dynamic structure factors. The presence of two 
distinct atomic species requires the calculation of the par- 
tial dynamic structure factors 

where summation over i and j run over all atoms of type 
a and /? respectively and Tj(t) is the position vector of 
atom j at time t. S a p(Q,E) is easily evaluated from 
the MD trajectories. In the NMA one introduces the 
displacements {u^} from the equilibrium positions {Xj} 
and expands these in terms of the normal modes; in the 
Q • Ui — ► limit the one excitation contribution can be 
calculated as 

S a p(Q,E) = k B Tj2 6(UJ ~ 2 UJX) x (1) 

A 

N a Np _ w 

VV ; = Q ■ e 4 (A)Q • c V :A^ x X 

where the Bose-Einstein factor has been approximated 
by its classical limit, Wij is the Debye- Waller factor and 
{lu\} and {ei(A)} are the eigenvalues and eigenvectors 
of the dynamical matrix of the system. Sn{Q,E) and 
Sx(Q,E) are then obtained by adding the partial dy- 
namic structure factors weighted by the atomic species 
concentrations and by the corresponding neutron or x- 
ray scattering lengths. An average over the independent 
directions of Q, i.e. Q = (±1, ±m, ±n)2it/L, with I, m, n 
integer numbers, is also performed. The reliability of the 



harmonic approximation is probed by the good agree- 
ment existing between the S(Q, E) calculated via NMA 
and via MD at finite temperature. Moreover, the systems 
with 648 and 5184 atonies give identical 5(Q, E) in their 
common Q range, indicating negligible size effects. 

Fig. 1 compares in a very satisfactory way INS data 
(squares) of ref. measured at 5.4 meV (corresponding 
to the position of the boson peak at T w 1100 K) with 
the neutron- weighted dynamic structure factor Sn(Q, E) 
calculated by our MD approach at the same temperature. 

We have also calculated the longitudinal current spec- 
tra C L {Q,E), defined as C L (Q,E) = E 2 S(Q,E)/Q 2 , at 
selected small Q values (Qmin — 1-47 nm _1 for the sy- 
stem with L = 42.78 A). These spectra, reported in Fig. 
2 for some low-Q values, have been fitted by the same 
model function (Damped Harmonic Oscillator, DHO) al- 
ready used J?J to analyse the experimental spectral shape. 
In Fig. 3 we report as full circles the best fit values of 
the excitation energies r2(Q), corresponding to the posi- 
tions of the current spectra maxima, and of the excita- 
tion widths r(Q). The form of the current spectra and 
the positions of the current spectra maxima are in qual- 
itative agreement with those of [^0|. The linear disper- 
sion observed for f2(Q) gives a longitudinal sound velocity 
which is slightly higher than that compatible with IXS 
data 0; this could be due to the thermal history and 
quench rate of the numerical model, which can produce 
a tempering of the vitreous structure. The evident linear 
Q dependence of the current spectra maxima positions, 
as already observable from the raw data of Fig. 2, in- 
dicates the propagating nature of the excitations up to 
w 20 meV, a range which includes the boson peak. 

This finding is not in contrast with earlier INS re- 
sults by Buchenau et al. pi] ], which demonstrate that 
pure plane-wave eigenvectors cannot explain the in- 
dependence of the S(Q, E) at E w 4 meV. Indeed, as 
it was shown also in other simulated glasses |22j, the 
eigenvectors of propagating modes present a random un- 
corrected component superimposed to the plane-wave. 

There is quantitative agreement between the Y(Q) va- 
lues obtained from our simulation and that measured by 
IXS for temperature ranging from 300 to 1800 K jL2|, 
both having a Y{Q) cx Q 2 behavior. The agreement of 
the r(Q) values found in our harmonic approximation 
with those measured in the real and simulated system, 
strongly indicates that the excitations linewidth, in the 
investigated Q range, is likely to be related to the struc- 
tural disorder rather than to anharmonicity or other dy- 
namical effects. As suggested in ref. R, the structural 
origin of the high frequency excitations widths is also 
supported by the temperature independence of T(Q) de- 
rived from IXS data. The picture becomes more compli- 
cated at small Q. Although the T(Q) ~ Q 2 law found 
at temperature above w 300 K by IXS and MD/NMA 
in the nm _1 range, quantitatively extrapolates over two 
decades in exchanged momentum down to the Brillouin 
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light scattering (BLS)range (Q w 0.02 nm _1 ), suggesting 
a common structural origin of the S(Q,u>) linewidth, at 
T < 100 K the T(Q) measured by BLS shows a strong 
temperature dependence p3[ . Consequently, a structural 
origin of T(Q) at high Q values is consistent with avail- 
able data, while at small Q the dynamical effects could 
be predominant. At present we have no explanation nei- 
ther for such different behaviour in different Q ranges 
nor for the T(Q) ~ Q 2 law observed at T >300 K in the 
Q = 0.02 5 imf 1 range. 

In ref. |l3fl , the localized nature of the high frequency 
excitations was derived from the lineshape analysis of 
the IXS data. However, the presence of a strong elas- 
tic contribution and the poor statistics, did not allow to 
discriminate between different lineshape models. A de- 
tailed lineshape analysis can be performed on the simu- 
lated Cl(Q,E) spectra. In the inset of Fig. 4 we report 
an example of fits of our calculated longitudinal current 
spectra with the expression proposed in ref. Re- 
liable fits can be obtained only leaving the important 
ujir parameter free. As reported in Fig. 4, the crossover 
frequency um, that according to the localization model 
should indicate the localization edge and therefore the 
boson peak energy, is strongly Q-dependent and assumes 
values definitely higher than that of the boson peak. If 
we impose the same value of the u>m parameter for the 
different Q spectra, the fits result unacceptably bad. 

To investigate the nature of the propagating excita- 
tions giving rise to the boson peak, we have also calcu- 
lated the transverse current spectra Ct(Q, E) at selected 
Q values. The transverse spectra are not experimentally 
accessible and are obtained from the MD via the corre- 
lation function of the transverse current or from NMA 
substituting (Q-) in Eq. |l| with (Qx). The transverse 
current spectra are also reported in Fig. 2, where we ob- 
serve the existence of a low energy excitation (E « 4 
meV at Q — 1.47 nm _1 ) which disperses with a sound 
velocity of w 4300 m/s up to « 15 mev at Q=5 nm _1 . 
The transverse modes become almost non-dispersing at 
E w 20 meV, an energy that is definitely higher than that 
characteristic of the boson peak. Also the propagating 
transverse dynamics may therefore contributes to the bo- 
son peak. It is worth to note that, in the Q range where 
the longitudinal and transverse modes show a propagat- 
ing character, the spectra Ct(Q,E) and Cl(Q,E) show 
a marked difference, again indicating the non-localized 
nature of the vibrational modes. Indeed, for localized 
modes one should not expect any difference between "lon- 
gitudinal" and "transverse" dynamics, being Q meaning- 
less. 

In conclusion, the numerical study of harmonic model 
of v-SiC>2 indicates that: i) The experimentally observed 
width of the excitations in the Q = 1.4-^3.5 nm _1 range 
is due to structural effects, i.e. to the fact that in this 
range Q is no longer a good quantum number in the 
topologically disordered silica. Anharmonicity or other 



dynamical effects (interaction with two level systems, soft 
potential modes, etc.) can only play a minor role, ii) 
Both longitudinal and transverse modes are propagating 
in the energy range covered by the boson peak: their 
dispersion relations saturate at much higher energy. Hi) 
The explanation of the numerical data with a cross-over 
model gives a Q dependent cross-over energy which is 
inconsistent with the model itself, iv) The transverse 
and longitudinal dynamics keep their differences up to 
at least Q — 5 nm" 1 , which indicates the non- localized 
nature of the vibrational states up to this Q-value. 

We greatefully akwnowledge F. Sette for many discus- 
sions and suggestions, P. Benassi and V. Mazzacurati for 
useful discussions, and U. Buchenau and S.R. Elliot for 
the preprints of ref. || and |l(J respectively. 
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FIGURE CAPTIONS 

Fig.l - Q dependence of the inelastic neutron scattering from vit- 
reous silica at 1.3 THz. The circle represent the data of 
ref. M at T = 1104K, the full line the results of the present 
NMA analysis multiplied by an arbitrary factor. 

Fig. 2 - X-ray longitudinal (line plus open circles) and transverse 
(line plus solid circles) current spectra obtained by the NMA 
at different Q values. Open and solid up-triangles refer to 
longitudinal and transverse spectra respectively, obtained by 
MD (216 Si0 2 units at 300 K). Full lines are the DHO fits 
of NMA longitudinal current spectra. 

Fig. 3 - a) Excitations energy Q(Q) from the DHO model for the 
IXS data at T = 1050K (open circles) and for the present 
NMA on the OK configuration, b) As m (a) but for the full 
width at half maximum of the excitations T(Q). 

Fig. 4 - Q-dependence of the cross over frequency loir resulting 
from the fits of the localised modes model of ref. |13] to 
the longitudinal current spectra, calculated from NMA. The 
inset shows an example of the fit quality. 
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